home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
The Datafile PD-CD 1 Issue 2
/
PDCD-1 - Issue 02.iso
/
_utilities
/
utilities
/
001
/
meschach
/
!Meschach
/
c
/
zfunc
< prev
next >
Wrap
Text File
|
1994-08-03
|
5KB
|
238 lines
/**************************************************************************
**
** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
**
** Meschach Library
**
** This Meschach Library is provided "as is" without any express
** or implied warranty of any kind with respect to this software.
** In particular the authors shall not be liable for any direct,
** indirect, special, incidental or consequential damages arising
** in any way from use of the software.
**
** Everyone is granted permission to copy, modify and redistribute this
** Meschach Library, provided:
** 1. All copies contain this copyright notice.
** 2. All modified copies shall carry a notice stating who
** made the last modification and the date of such modification.
** 3. No charge is made for this software or works derived from it.
** This clause shall not be construed as constraining other software
** distributed on the same medium as this software, nor is a
** distribution fee considered a charge.
**
***************************************************************************/
/*
Elementary functions for complex numbers
-- if not already defined
*/
#include <math.h>
#include "zmatrix.h"
static char rcsid[] = "$Id: zfunc.c,v 1.1 1994/01/13 04:28:29 des Exp $";
#ifndef COMPLEX_H
#ifndef zmake
/* zmake -- create complex number real + i*imag */
complex zmake(real,imag)
double real, imag;
{
complex w; /* == real + i*imag */
w.re = real; w.im = imag;
return w;
}
#endif
#ifndef zneg
/* zneg -- returns negative of z */
complex zneg(z)
complex z;
{
z.re = - z.re;
z.im = - z.im;
return z;
}
#endif
#ifndef zabs
/* zabs -- returns |z| */
double zabs(z)
complex z;
{
Real x, y, tmp;
int x_expt, y_expt;
/* Note: we must ensure that overflow does not occur! */
x = ( z.re >= 0.0 ) ? z.re : -z.re;
y = ( z.im >= 0.0 ) ? z.im : -z.im;
if ( x < y )
{
tmp = x;
x = y;
y = tmp;
}
x = frexp(x,&x_expt);
y = frexp(y,&y_expt);
y = ldexp(y,y_expt-x_expt);
tmp = sqrt(x*x+y*y);
return ldexp(tmp,x_expt);
}
#endif
#ifndef zadd
/* zadd -- returns z1+z2 */
complex zadd(z1,z2)
complex z1, z2;
{
complex z;
z.re = z1.re + z2.re;
z.im = z1.im + z2.im;
return z;
}
#endif
#ifndef zsub
/* zsub -- returns z1-z2 */
complex zsub(z1,z2)
complex z1, z2;
{
complex z;
z.re = z1.re - z2.re;
z.im = z1.im - z2.im;
return z;
}
#endif
#ifndef zmlt
/* zmlt -- returns z1*z2 */
complex zmlt(z1,z2)
complex z1, z2;
{
complex z;
z.re = z1.re * z2.re - z1.im * z2.im;
z.im = z1.re * z2.im + z1.im * z2.re;
return z;
}
#endif
#ifndef zinv
/* zmlt -- returns 1/z */
complex zinv(z)
complex z;
{
Real x, y, tmp;
int x_expt, y_expt;
if ( z.re == 0.0 && z.im == 0.0 )
error(E_SING,"zinv");
/* Note: we must ensure that overflow does not occur! */
x = ( z.re >= 0.0 ) ? z.re : -z.re;
y = ( z.im >= 0.0 ) ? z.im : -z.im;
if ( x < y )
{
tmp = x;
x = y;
y = tmp;
}
x = frexp(x,&x_expt);
y = frexp(y,&y_expt);
y = ldexp(y,y_expt-x_expt);
/* 'ldexp' call moved to 'tmp' assignment by Paul Field */
tmp = (1.0/(x*x + y*y)) * ldexp(1.0,-2*x_expt);
z.re = z.re*tmp /* * ldexp(1.0,-2*x_expt) */;
z.im = -z.im*tmp /* * ldexp(1.0,-2*x_expt) */;
return z;
}
#endif
#ifndef zdiv
/* zdiv -- returns z1/z2 */
complex zdiv(z1,z2)
complex z1, z2;
{
return zmlt(z1,zinv(z2));
}
#endif
#ifndef zsqrt
/* zsqrt -- returns sqrt(z); uses branch with Re sqrt(z) >= 0 */
complex zsqrt(z)
complex z;
{
complex w; /* == sqrt(z) at end */
Real alpha;
alpha = sqrt(0.5*(fabs(z.re) + zabs(z)));
if ( z.re >= 0.0 )
{
w.re = alpha;
w.im = z.im / (2.0*alpha);
}
else
{
w.re = fabs(z.im)/(2.0*alpha);
w.im = ( z.im >= 0 ) ? alpha : - alpha;
}
return w;
}
#endif
#ifndef zexp
/* zexp -- returns exp(z) */
complex zexp(z)
complex z;
{
complex w; /* == exp(z) at end */
Real r;
r = exp(z.re);
w.re = r*cos(z.im);
w.im = r*sin(z.im);
return w;
}
#endif
#ifndef zlog
/* zlog -- returns log(z); uses principal branch with -pi <= Im log(z) <= pi */
complex zlog(z)
complex z;
{
complex w; /* == log(z) at end */
w.re = log(zabs(z));
w.im = atan2(z.im,z.re);
return w;
}
#endif
#ifndef zconj
complex zconj(z)
complex z;
{
complex w; /* == conj(z) */
w.re = z.re;
w.im = - z.im;
return w;
}
#endif
#endif